Let’s import at all of the data that Simon and Matt pulled:

library(tidyverse)
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
── Attaching packages ───────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.1     ✔ purrr   0.3.2
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   1.0.0     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine
library(maps)

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map
library(ggthemes)


DeltasClean <- read_csv("../data/out/deltas_clean_v2.csv") 
Parsed with column specification:
cols(
  Delta = col_character(),
  location = col_character(),
  surface = col_character(),
  year = col_double(),
  month = col_double(),
  ndvi = col_double(),
  red = col_double(),
  evi = col_double(),
  savi = col_double(),
  gr = col_double(),
  ndssi = col_double()
)
DeltaLocations <- read_csv("../data/DeltaLocations.csv")
Parsed with column specification:
cols(
  Deltas = col_character(),
  Lat = col_double(),
  Lon = col_double()
)

As a reminder, for each of the 47 deltas there are measurements of Land & Water areas at Upstream, Downstream and ‘Middle’ locations on the delta. We first lump all the observations together, and look to see which Deltas have many observations:

#counts per delta
DeltaCounts <- count(DeltasClean, Delta)
DeltaCounts

Now, by each month.. where the colorbar represents the number of observations (n) for each month for a given delta:

DeltaObsPerMonth <- count(DeltasClean, Delta, month)

ggplot(DeltaObsPerMonth, aes(y = Delta, x = month, fill=n)) + geom_tile() + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
  expand_limits(x = c(1,12)) + 
  scale_fill_gradient( trans = 'log' )

In the above heat map, dark colors (and no color) represent data paucity (and data gaps). Deltas with light colors (e.g., the Parana, Nile, Ebro, Colorado, Brahmani) have lots of data, spread out through the months of the year.

I’ll remove/subset the deltas with sparse coverage (specifically, months with no coverage)….


# need 10 data points per month for NDSSI and NDVI
EnoughObsPerMonth <- DeltasClean %>% ungroup() %>%
  count(Delta, month, surface) %>% 
  group_by(surface) %>%
  filter( n >= 5)

#find deltas missing a given month of observations
DeltaMonthCounts <- EnoughObsPerMonth %>%
  ungroup() %>%
  count(Delta)

# need 12 months of water and land obs, so 24 mo total
EnoughMonths <- DeltaMonthCounts %>%
 filter( n == 24)

CompleteObsDeltas <- pull(EnoughMonths, Delta)

#remove them
DeltasCleaner <- DeltasClean %>%
  filter(Delta %in% CompleteObsDeltas)

#add the real dates in month date format
DeltasCleaner$date <- as.Date(paste(DeltasCleaner$year, DeltasCleaner$month, "01", sep="-"), "%Y-%m-%d")

#remove intermediate data
rm(CompleteObsDeltas, EnoughMonths, DeltaMonthCounts)

#EnoughMonths

and extract some metrics; specifically I will make a timeseries of NDVI and NDSSI for each delta using the mean value for each month.

From a quick look at the dat (not shown until much later in the doc):

#take the mean NDVI and NDSSI for each month, for each delta
DeltaMeans <- DeltasCleaner %>%
  group_by(Delta, month, surface) %>%
  summarize(MeanNDVI = mean(ndvi, na.rm = TRUE), MeanNDSSI = mean(ndssi, na.rm = TRUE))

#make a 9 column data frame with:
#delta, 
#max and min NDVI month, 
#NDSSI max and min month, 
#max and min values for both NDVI and NDSSI

#####

DeltaMaxNDVI <- 
  DeltaMeans %>% 
  filter(surface == 'Land')  %>% 
  select (-c(MeanNDSSI)) %>% 
  group_by(Delta) %>% 
  slice(which.max(MeanNDVI)) %>% 
  rename(MaxMeanNDVImonth = month, MaxMeanNDVI = MeanNDVI)

DeltaMaxNDSSI <- 
  DeltaMeans %>% 
  filter(surface == 'Water')  %>% 
  select (-c(MeanNDVI)) %>% 
  group_by(Delta) %>% 
  slice(which.max(MeanNDSSI)) %>% 
  rename(MaxMeanNDSSImonth = month, MaxMeanNDSSI = MeanNDSSI)


DeltaMinNDVI <- 
  DeltaMeans %>% 
  filter(surface == 'Land')  %>% 
  select (-c(MeanNDSSI)) %>% 
  group_by(Delta) %>% 
  slice(which.min(MeanNDVI)) %>% 
  rename(MinMeanNDVImonth = month, MinMeanNDVI = MeanNDVI)

DeltaMinNDSSI <- 
  DeltaMeans %>% 
  filter(surface == 'Water')  %>% 
  select (-c(MeanNDVI)) %>% 
  group_by(Delta) %>% 
  slice(which.min(MeanNDSSI)) %>% 
  rename(MinMeanNDSSImonth = month, MinMeanNDSSI = MeanNDSSI)


#join into 1 dataframe
DeltaMaxMin <- left_join(DeltaMaxNDVI, DeltaMaxNDSSI, by = 'Delta') %>% 
  left_join(.,DeltaMinNDVI, by = 'Delta') %>% 
  left_join(.,DeltaMinNDSSI, by = 'Delta')

#remove intermediate data
rm(DeltaMaxNDVI, DeltaMaxNDSSI, DeltaMinNDSSI,DeltaMinNDVI)

And now we can look at the phase shifts between these two time series (the timeseries of mean NDVI and mean NDSSI). Here are the phase shifts (in month) for each delta:

#compare offset
DeltaMaxMin <- mutate(DeltaMaxMin, 
                      MinOffset = if_else(abs(MinMeanNDVImonth - MinMeanNDSSImonth) > 6, 
                                          12 - abs(MinMeanNDVImonth - MinMeanNDSSImonth),
                                          abs(MinMeanNDVImonth - MinMeanNDSSImonth)),
                      MaxOffset = if_else(abs(MaxMeanNDVImonth - MaxMeanNDSSImonth) > 6, 
                                         12 - abs(MaxMeanNDVImonth - MaxMeanNDSSImonth),
                                          abs(MaxMeanNDVImonth - MaxMeanNDSSImonth)),
                      OffsetDiff = abs(MaxOffset - MinOffset),
                      rangeNDVI = (MaxMeanNDVI - MinMeanNDVI), 
                      rangeNDSSI = (MaxMeanNDSSI - MinMeanNDSSI),
                      halfPeriodNDVI = if_else(abs(MaxMeanNDVImonth - MinMeanNDVImonth) > 6, 
                                          12 - abs(MaxMeanNDVImonth - MinMeanNDVImonth),
                                          abs(MaxMeanNDVImonth - MinMeanNDVImonth)),
                      halfPeriodNDSSI = if_else(abs(MaxMeanNDSSImonth - MinMeanNDSSImonth) > 6, 
                                          12 - abs(MaxMeanNDSSImonth - MinMeanNDSSImonth),
                                          abs(MaxMeanNDSSImonth - MinMeanNDSSImonth)), )

# DeltaMaxMin <- 
#   DeltaMaxMin   %>%
#   select (c(Delta, MinOffset, MaxOffset, OffsetDiff, rangeNDVI, rangeNDSSI,MaxMeanNDSSI,MinMeanNDSSI,MaxMeanNDVI,MinMeanNDVI)) 

DeltaMaxMin

ggplot(DeltaMaxMin, aes(y = Delta, x = MaxOffset)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("MaxOffset")


ggplot(DeltaMaxMin, aes(y = Delta, x = MinOffset)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("MinOffset")


ggplot(DeltaMaxMin, aes(y = Delta, x = OffsetDiff)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("Offset Difference")

ggplot(DeltaMaxMin, aes(x = MaxMeanNDVImonth)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.5) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) + 
  labs(x = "Month") +
  ggtitle("Month of maximum mean NDVI")


ggplot(DeltaMaxMin, aes(x = MaxMeanNDSSImonth)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.5) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) + 
  labs(x = "Month") +
  ggtitle("Month of maximum mean NDSSI")


ggplot(DeltaMaxMin, aes(x = MaxOffset)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("Months Offset between NDVI and NDSSI")


ggplot(DeltaMaxMin, aes(x = halfPeriodNDVI)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("half period NDVI ")


ggplot(DeltaMaxMin, aes(x = halfPeriodNDSSI)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("half period NDSSI")

#extract NDVI value for each delta a the month of max NDSSI value

DeltaNDVIatMaxNDSSI <- DeltaMaxMin %>%
  select(Delta,MaxMeanNDSSImonth)

DeltaMeansToJoin <- DeltaMeans %>%
  filter(surface == 'Land')

DeltaNDVIatMaxNDSSI <- left_join(DeltaNDVIatMaxNDSSI, DeltaMeansToJoin, 
                                 by = c('Delta', 'MaxMeanNDSSImonth' ='month'))
 
DeltaNDVIatMaxNDSSI <- DeltaNDVIatMaxNDSSI %>%
  select (-c(surface, MeanNDSSI))

DeltaNDVIatMaxNDSSI

ggplot(DeltaNDVIatMaxNDSSI, aes(x = MeanNDVI)) + 
  geom_dotplot(binwidth = 0.1, dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL)  + 
  xlim(0,1) +
  labs(x = "NDVI") +
  ggtitle("NDVI at month of maximum mean NDSSI")

#extract NDVI value for each delta a the month of min NDSSI value

DeltaNDVIatMinNDSSI <- DeltaMaxMin %>%
  select(Delta,MinMeanNDSSImonth)

DeltaMeansToJoin <- DeltaMeans %>%
  filter(surface == 'Land')

DeltaNDVIatMinNDSSI <- left_join(DeltaNDVIatMinNDSSI, DeltaMeansToJoin, 
                                 by = c('Delta', 'MinMeanNDSSImonth' ='month'))
 
DeltaNDVIatMinNDSSI <- DeltaNDVIatMinNDSSI %>%
  select (-c(surface, MeanNDSSI))

DeltaNDVIatMaxNDSSI

ggplot(DeltaNDVIatMinNDSSI, aes(x = MeanNDVI)) + 
  geom_dotplot(binwidth = 0.1, dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL)  + 
  xlim(0,1) +
  labs(x = "NDVI") +
  ggtitle("NDVI at month of min mean NDSSI")

And here are the phase shifts/offsets against other measured parameters for each delta. The range, max and mean of NDVI and NDSSI is calculated from the timeseries, so it is really the max, min, and range of the monthly means (i.e., the maximum of the means, the minimum of the means, and the range of the mean). Offset is measured in months.

slice1 <- ggplot(DeltaMaxMin, aes(y = rangeNDVI, x = MaxOffset)) + geom_point() 
# + scale_x_discrete(limits = c(1:6), breaks = c(1:6)) + expand_limits(x = c(1,6)) 

slice2 <- ggplot(DeltaMaxMin, aes(y = rangeNDVI, x = rangeNDSSI)) + geom_point() 
slice3 <- ggplot(DeltaMaxMin, aes(y = rangeNDSSI, x = MaxOffset)) + geom_point() 

slice4 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = rangeNDVI)) + geom_point() 
slice5 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = rangeNDSSI)) + geom_point() 
slice6 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = MaxOffset)) + geom_point() 

slice7 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = MaxMeanNDVI)) + geom_point() 
slice8 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = rangeNDSSI)) + geom_point() 
slice9 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = MaxOffset)) + geom_point() 
slice10 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = rangeNDVI)) + geom_point() 

grid.arrange(slice1, slice2, slice3, slice4, slice5, slice6, ncol=3)

grid.arrange(slice7, slice8, slice9, slice10, ncol=2)


#remove those panels
rm(slice1, slice2, slice3, slice4, slice5, slice6,slice7, slice8, slice9, slice10)

Join Latitude and longitude data

DeltaDatawLocations <- left_join(DeltaMaxMin, DeltaLocations, by = c("Delta" = "Deltas"))

DeltaDatawLocations <- DeltaDatawLocations %>%
  mutate(Absolute_Latitude= abs(Lat))

plot params vs lat

loc1 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxOffset)) + geom_point() 
loc2 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = rangeNDSSI)) + geom_point() 
loc3 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = rangeNDVI)) + geom_point() 
loc4 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxMeanNDVI)) + geom_point() 
loc5 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxMeanNDSSI)) + geom_point() 

grid.arrange(loc1, loc2, loc3, loc4, loc5, ncol=2)


loc1

#ggsave("loc1.pdf", width = 4, height = 4)
loc2

#ggsave("loc2.pdf", width = 4, height = 4)
loc3

#ggsave("loc3.pdf", width = 4, height = 4)

#remove those panels
rm(loc1, loc2, loc3, loc4, loc5)
#find the linear model 
DeltaOffset_lm <- lm( Absolute_Latitude ~ MaxOffset, data = DeltaDatawLocations) 

summary(DeltaOffset_lm)

Call:
lm(formula = Absolute_Latitude ~ MaxOffset, data = DeltaDatawLocations)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.779  -8.983   0.775   8.377  20.772 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   10.820      4.629   2.337  0.02653 * 
MaxOffset      4.079      1.157   3.525  0.00143 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.06 on 29 degrees of freedom
Multiple R-squared:  0.2999,    Adjusted R-squared:  0.2758 
F-statistic: 12.43 on 1 and 29 DF,  p-value: 0.001428
ggplot(DeltaDatawLocations, aes(x = Absolute_Latitude, y = MaxOffset)) + 
  geom_point() +
  geom_smooth(mapping = aes(x = Absolute_Latitude, y = MaxOffset, ), method=lm ) 

NA
NA
NA

Now for the maps:

world <- ggplot() +
  borders("world", colour = "gray85", fill = "gray80") +
  theme_map() 

DeltaOffsetMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = MaxOffset),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") +
  ggtitle("Offset Between NDVI peak on Land and NDSSI peak in water")

DeltaOffsetMap

#ggsave("DeltaOffsetMap.pdf", width = 6, height = 4)
world <- ggplot() +
  borders("world", colour = "gray85", fill = "gray80") +
  theme_map() 

DeltaNDVIrangeMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = rangeNDVI),
             data = DeltaDatawLocations,
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") + 
  ggtitle("NDVI range")


DeltaNDSSIrangeMap  <- world +
  geom_point(aes(x = Lon, y = Lat, color = rangeNDSSI),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low = "yellow") + ggtitle("NDSSI range") 

DeltaNDVIrangeMap 

#ggsave("DeltaNDVIrangeMap.pdf", width = 6, height = 4)
DeltaNDSSIrangeMap

#ggsave("DeltaNDSSIrangeMap.pdf", width = 6, height = 4)

Let’s look at some of the data. To quantify the water, we use NDSSI. to quantify land, we use NDVI.here is the time series for the Parana.

DeltaPlotter <- function(DeltaName) {
  #Counts each month
  numVeg <- DeltasCleaner %>%
    select(Delta, surface, month, ndvi) %>%
    filter(Delta == DeltaName & surface == "Land" & !is.na(ndvi)) %>%
    group_by(month) %>%
    summarize(n = n())
  
  numSed <- DeltasCleaner %>%
    select(Delta, surface, month, ndssi) %>%
    filter(Delta == DeltaName &
             surface == "Water" & !is.na(ndssi)) %>%
    group_by(month) %>%
    summarize(n = n())
  
  #Highlight the Maximum and Minimum Month for each delta, NDVI and NDSSI
  
  #LAND
  Veg <-
    ggplot(data = filter(DeltasCleaner, Delta == DeltaName &
                           surface == "Land")) +
    geom_boxplot(aes(x = month, y = ndvi, group = month)) +
    scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
    expand_limits(x = c(1, 12)) +
    ggtitle(DeltaName) +
    #geom_text(data = numVeg, aes(y = 1.05, x = month, label = n)) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Land" & month == DeltaMaxMin$MaxMeanNDVImonth[DeltaMaxMin$Delta == DeltaName] 
      ),
      aes(x = month, y = ndvi, group = month),
      fill = "green"
    ) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName & 
          surface == "Land" & month == DeltaMaxMin$MinMeanNDVImonth[DeltaMaxMin$Delta ==DeltaName]
      ),
      aes(x = month, y = ndvi, group = month),
      fill = "blue"
    )
  
  
  Sed <-
    ggplot(data = filter(DeltasCleaner, Delta == DeltaName &
                           surface == "Water")) +
    geom_boxplot(aes(x = month, y = ndssi, group = month)) +
    scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
    expand_limits(x = c(1, 12)) +
    #geom_text(data = numSed, aes(y = 1.05, x = month, label = n)) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Water" & month == DeltaMaxMin$MaxMeanNDSSImonth[DeltaMaxMin$Delta == DeltaName]
      ),
      aes(x = month, y = ndssi, group = month),
      fill = "green"
    ) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Water" & month == DeltaMaxMin$MinMeanNDSSImonth[DeltaMaxMin$Delta == DeltaName]
      ),
      aes(x = month, y = ndssi, group = month),
      fill = "blue"
    )
  
  return(grid.arrange(Veg, Sed, nrow = 2))
}

Here is an ‘Offset’ example from the Parana:

DeltaPlotter("Parana")

and now an in-phase example from

DeltaPlotter("Magdalena")

and the Ebro:

DeltaPlotter("Ebro")

and the Tana

DeltaPlotter("Tana")

and the Tana

DeltaPlotter("Godavari")

DeltaPlotter("Krishna")

Now for some work with GRDC data:

#import the data (monthly means for 21 stations)
DeltasGRDC  <- read_csv("../data/GRDCstations.csv")
Parsed with column specification:
cols(
  Deltas = col_character(),
  GRDC_Station = col_double(),
  Time_Series_Length = col_character(),
  January = col_double(),
  February = col_double(),
  March = col_double(),
  April = col_double(),
  May = col_double(),
  June = col_double(),
  July = col_double(),
  August = col_double(),
  September = col_double(),
  October = col_double(),
  November = col_double(),
  December = col_double()
)
#calculate the mean of the monthly means
DeltasGRDCwMean <- DeltasGRDC %>% 
    rowwise() %>% 
    mutate(MMD=mean(c(January,February,March,April,
                    May,June,July,August,
                    September,October,November,December)))

DeltasGRDCwMean <- DeltasGRDCwMean %>% 
  rowwise() %>% 
  mutate(Range_Discharge = max(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December)) - 
           min(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December))
           )
        

#join to location data:
DeltawLocGRDC <- left_join(DeltaDatawLocations, DeltasGRDCwMean, by = c("Delta" = "Deltas"))


#plot mean of monthly means against NDSSI
ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = rangeNDSSI)) + geom_point() + scale_y_continuous(trans='log10')


#ggsave("GRDCNDSSI.pdf", width = 6, height = 4)

ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = MaxMeanNDSSI)) + geom_point() + scale_y_continuous(trans='log10')

#import the data (monthly means for 21 stations)
DeltasGRDC  <- read_csv("../data/GRDCstations.csv")
Parsed with column specification:
cols(
  Deltas = col_character(),
  GRDC_Station = col_double(),
  Time_Series_Length = col_character(),
  January = col_double(),
  February = col_double(),
  March = col_double(),
  April = col_double(),
  May = col_double(),
  June = col_double(),
  July = col_double(),
  August = col_double(),
  September = col_double(),
  October = col_double(),
  November = col_double(),
  December = col_double()
)
#calculate the mean of the monthly means
DeltasGRDCwMean <- DeltasGRDC %>% 
    rowwise() %>% 
    mutate(MMD=mean(c(January,February,March,April,
                    May,June,July,August,
                    September,October,November,December)))

DeltasGRDCwMean <- DeltasGRDCwMean %>% 
  rowwise() %>% 
  mutate(Range_Discharge = max(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December)) - 
           min(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December))
           )
        

#join to location data:
DeltawLocGRDC <- left_join(DeltaDatawLocations, DeltasGRDCwMean, by = c("Delta" = "Deltas"))


#plot mean of monthly means against NDSSI
ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = rangeNDSSI)) + geom_point() + scale_y_continuous(trans='log10')


#ggsave("GRDCNDSSI.pdf", width = 6, height = 4)
ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = MaxMeanNDSSI)) + geom_point() + scale_y_continuous(trans='log10')

#rename the months by numbers and tidy the GRDC data
DeltasDischarge <- DeltasGRDC %>%
  rename(Delta = Deltas,"1" = January, "2"= February, "3"= March, "4"= April,
         "5"=May, "6"=June, "7"=July, "8"= August, "9" = September, "10"=October, 
         "11"=November, "12"=December) %>%
  select(Delta, "1" , "2" , "3", "4","5", "6", "7", "8", "9", "10", "11", "12") %>%
  pivot_longer(-Delta, names_to = "month", values_to = "discharge")

#find max GRDC month for each delta
DeltaMaxDischarge <- 
  DeltasDischarge %>% 
  group_by(Delta) %>% 
  slice(which.max(discharge)) %>% 
  rename(MaxDischargeMonth = month, MaxDischarge = discharge) 


DeltaMaxDischarge$MaxDischargeMonth = as.numeric(DeltaMaxDischarge$MaxDischargeMonth)


#join with other delta data
DeltaMaxMinDischarge <- left_join(DeltaMaxMin, DeltaMaxDischarge, by = 'Delta')

#calculate offset
DeltaMaxMinDischarge <- DeltaMaxMinDischarge %>%
  mutate( DissOff = if_else(abs(MaxDischargeMonth - MaxMeanNDSSImonth) > 6,
                            12 - abs(MaxDischargeMonth - MaxMeanNDSSImonth),
                            abs(MaxDischargeMonth - MaxMeanNDSSImonth))
          )

#Compare offset with NDSSI (deltamaxmin$maxmeanNDSSImonth)
ggplot(DeltaMaxMinDischarge, aes(y = Delta, x = DissOff)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("DisOffset")

#join lat data
DeltaDatawLocations <- left_join(DeltaMaxMinDischarge, DeltaLocations, by = c("Delta" = "Deltas"))

DeltaDatawLocations <- DeltaDatawLocations %>%
  mutate(Absolute_Latitude= abs(Lat))

# plot offset on graph by lat
ggplot(DeltaDatawLocations, aes(x = Absolute_Latitude, y = DissOff)) + geom_point() +
  scale_color_gradient(low = "yellow", high = "red", na.value = NA) 

  #+ geom_smooth(mapping = aes(x = Absolute_Latitude, y = DissOff, ), method=lm ) 

#ggsave("DisOffset.pdf", width = 6, height = 4)
#plot offset on map
DeltaDisOffsetMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = DissOff),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") +
  ggtitle("Offset Between GRDC discharge peak and NDSSI peak in water")

DeltaDisOffsetMap

#ggsave("DeltaOffsetMap.pdf", width = 6, height = 4)
LS0tCnRpdGxlOiAiRGVsdGFzIE5vdGVib29rIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIGRmX3ByaW50OiBwYWdlZAogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKICB3b3JkX2RvY3VtZW50OiBkZWZhdWx0CmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKTGV0J3MgaW1wb3J0IGF0IGFsbCBvZiB0aGUgZGF0YSB0aGF0IFNpbW9uIGFuZCBNYXR0IHB1bGxlZDoKCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShncmlkRXh0cmEpCmxpYnJhcnkobWFwcykKbGlicmFyeShnZ3RoZW1lcykKCgpEZWx0YXNDbGVhbiA8LSByZWFkX2NzdigiLi4vZGF0YS9vdXQvZGVsdGFzX2NsZWFuX3YyLmNzdiIpIApEZWx0YUxvY2F0aW9ucyA8LSByZWFkX2NzdigiLi4vZGF0YS9EZWx0YUxvY2F0aW9ucy5jc3YiKQpgYGAKCkFzIGEgcmVtaW5kZXIsIGZvciBlYWNoIG9mIHRoZSA0NyBkZWx0YXMgdGhlcmUgYXJlIG1lYXN1cmVtZW50cyBvZiBMYW5kICYgV2F0ZXIgYXJlYXMgYXQgVXBzdHJlYW0sIERvd25zdHJlYW0gYW5kICdNaWRkbGUnIGxvY2F0aW9ucyBvbiB0aGUgZGVsdGEuIFdlIGZpcnN0IGx1bXAgYWxsIHRoZSBvYnNlcnZhdGlvbnMgdG9nZXRoZXIsIGFuZCBsb29rIHRvIHNlZSB3aGljaCBEZWx0YXMgaGF2ZSBtYW55IG9ic2VydmF0aW9uczoKCmBgYHtyfQojY291bnRzIHBlciBkZWx0YQpEZWx0YUNvdW50cyA8LSBjb3VudChEZWx0YXNDbGVhbiwgRGVsdGEpCkRlbHRhQ291bnRzCmBgYAoKCk5vdywgYnkgZWFjaCBtb250aC4uIHdoZXJlIHRoZSBjb2xvcmJhciByZXByZXNlbnRzIHRoZSBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zIChuKSBmb3IgZWFjaCBtb250aCBmb3IgYSBnaXZlbiBkZWx0YToKYGBge3J9CkRlbHRhT2JzUGVyTW9udGggPC0gY291bnQoRGVsdGFzQ2xlYW4sIERlbHRhLCBtb250aCkKCmdncGxvdChEZWx0YU9ic1Blck1vbnRoLCBhZXMoeSA9IERlbHRhLCB4ID0gbW9udGgsIGZpbGw9bikpICsgZ2VvbV90aWxlKCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMToxMiksIGJyZWFrcyA9IGMoMToxMikpICsKICBleHBhbmRfbGltaXRzKHggPSBjKDEsMTIpKSArIAogIHNjYWxlX2ZpbGxfZ3JhZGllbnQoIHRyYW5zID0gJ2xvZycgKQoKYGBgCgoKSW4gdGhlIGFib3ZlIGhlYXQgbWFwLCBkYXJrIGNvbG9ycyAoYW5kIG5vIGNvbG9yKSByZXByZXNlbnQgZGF0YSBwYXVjaXR5IChhbmQgZGF0YSBnYXBzKS4gRGVsdGFzIHdpdGggbGlnaHQgY29sb3JzIChlLmcuLCB0aGUgUGFyYW5hLCBOaWxlLCBFYnJvLCBDb2xvcmFkbywgQnJhaG1hbmkpIGhhdmUgbG90cyBvZiBkYXRhLCBzcHJlYWQgb3V0IHRocm91Z2ggdGhlIG1vbnRocyBvZiB0aGUgeWVhci4KCgpJJ2xsIHJlbW92ZS9zdWJzZXQgdGhlIGRlbHRhcyB3aXRoIHNwYXJzZSBjb3ZlcmFnZSAoc3BlY2lmaWNhbGx5LCBtb250aHMgd2l0aCBubyBjb3ZlcmFnZSkuLi4uIApgYGB7cn0KCiMgbmVlZCAxMCBkYXRhIHBvaW50cyBwZXIgbW9udGggZm9yIE5EU1NJIGFuZCBORFZJCkVub3VnaE9ic1Blck1vbnRoIDwtIERlbHRhc0NsZWFuICU+JSB1bmdyb3VwKCkgJT4lCiAgY291bnQoRGVsdGEsIG1vbnRoLCBzdXJmYWNlKSAlPiUgCiAgZ3JvdXBfYnkoc3VyZmFjZSkgJT4lCiAgZmlsdGVyKCBuID49IDUpCgojZmluZCBkZWx0YXMgbWlzc2luZyBhIGdpdmVuIG1vbnRoIG9mIG9ic2VydmF0aW9ucwpEZWx0YU1vbnRoQ291bnRzIDwtIEVub3VnaE9ic1Blck1vbnRoICU+JQogIHVuZ3JvdXAoKSAlPiUKICBjb3VudChEZWx0YSkKCiMgbmVlZCAxMiBtb250aHMgb2Ygd2F0ZXIgYW5kIGxhbmQgb2JzLCBzbyAyNCBtbyB0b3RhbApFbm91Z2hNb250aHMgPC0gRGVsdGFNb250aENvdW50cyAlPiUKIGZpbHRlciggbiA9PSAyNCkKCkNvbXBsZXRlT2JzRGVsdGFzIDwtIHB1bGwoRW5vdWdoTW9udGhzLCBEZWx0YSkKCiNyZW1vdmUgdGhlbQpEZWx0YXNDbGVhbmVyIDwtIERlbHRhc0NsZWFuICU+JQogIGZpbHRlcihEZWx0YSAlaW4lIENvbXBsZXRlT2JzRGVsdGFzKQoKI2FkZCB0aGUgcmVhbCBkYXRlcyBpbiBtb250aCBkYXRlIGZvcm1hdApEZWx0YXNDbGVhbmVyJGRhdGUgPC0gYXMuRGF0ZShwYXN0ZShEZWx0YXNDbGVhbmVyJHllYXIsIERlbHRhc0NsZWFuZXIkbW9udGgsICIwMSIsIHNlcD0iLSIpLCAiJVktJW0tJWQiKQoKI3JlbW92ZSBpbnRlcm1lZGlhdGUgZGF0YQpybShDb21wbGV0ZU9ic0RlbHRhcywgRW5vdWdoTW9udGhzLCBEZWx0YU1vbnRoQ291bnRzKQoKI0Vub3VnaE1vbnRocwpgYGAKCmFuZCBleHRyYWN0IHNvbWUgbWV0cmljczsgc3BlY2lmaWNhbGx5IEkgd2lsbCBtYWtlIGEgdGltZXNlcmllcyBvZiBORFZJIGFuZCBORFNTSSBmb3IgZWFjaCBkZWx0YSB1c2luZyB0aGUgbWVhbiB2YWx1ZSBmb3IgZWFjaCBtb250aC4KCkZyb20gYSBxdWljayBsb29rIGF0IHRoZSBkYXQgKG5vdCBzaG93biB1bnRpbCBtdWNoIGxhdGVyIGluIHRoZSBkb2MpOgoKKiBQZWFrcyBpbiBORFZJIGFyZSBub3QgYWx3YXlzIGNvcnJlbGF0ZWQgd2l0aCBwZWFrcyBpbiBORFNTSS4gaW4gb3RoZXIgd29yZHMsIG1heGltdW1zIGluIHNlZCBjb25jZW50cmF0aW9uIGFyZSBub3QgZHVyaW5nIG1heGltdW1zIGluIHN0YW5kaW5nIGJpb21hc3MuCiogVGhlIHBlYWtzIGluIGJvdGggdGltZXNlcmllcyBzaGlmdCBhcm91bmQgZGVwZW5kaW5nIG9uIHRoZSBkZWx0YToKICAgKyBsb29rIGF0IHRoZSBjb3JyZWxhdGlvbiBpbiB0aGUgTmlsZS4KICAgKyBUaGUgYW50aWNvcnJlbGF0aW9uIGluIHRoZSBFYnJvLAogICArIFRoZSBzbGlnaHQgcGhhc2Ugc2hpZnQgaW4gdGhlIFBhcmFuYS4KCgpgYGB7ciAgaW5jbHVkZSA9IFRSVUV9CiN0YWtlIHRoZSBtZWFuIE5EVkkgYW5kIE5EU1NJIGZvciBlYWNoIG1vbnRoLCBmb3IgZWFjaCBkZWx0YQpEZWx0YU1lYW5zIDwtIERlbHRhc0NsZWFuZXIgJT4lCiAgZ3JvdXBfYnkoRGVsdGEsIG1vbnRoLCBzdXJmYWNlKSAlPiUKICBzdW1tYXJpemUoTWVhbk5EVkkgPSBtZWFuKG5kdmksIG5hLnJtID0gVFJVRSksIE1lYW5ORFNTSSA9IG1lYW4obmRzc2ksIG5hLnJtID0gVFJVRSkpCgojbWFrZSBhIDkgY29sdW1uIGRhdGEgZnJhbWUgd2l0aDoKI2RlbHRhLCAKI21heCBhbmQgbWluIE5EVkkgbW9udGgsIAojTkRTU0kgbWF4IGFuZCBtaW4gbW9udGgsIAojbWF4IGFuZCBtaW4gdmFsdWVzIGZvciBib3RoIE5EVkkgYW5kIE5EU1NJCgojIyMjIwoKRGVsdGFNYXhORFZJIDwtIAogIERlbHRhTWVhbnMgJT4lIAogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykgICU+JSAKICBzZWxlY3QgKC1jKE1lYW5ORFNTSSkpICU+JSAKICBncm91cF9ieShEZWx0YSkgJT4lIAogIHNsaWNlKHdoaWNoLm1heChNZWFuTkRWSSkpICU+JSAKICByZW5hbWUoTWF4TWVhbk5EVkltb250aCA9IG1vbnRoLCBNYXhNZWFuTkRWSSA9IE1lYW5ORFZJKQoKRGVsdGFNYXhORFNTSSA8LSAKICBEZWx0YU1lYW5zICU+JSAKICBmaWx0ZXIoc3VyZmFjZSA9PSAnV2F0ZXInKSAgJT4lIAogIHNlbGVjdCAoLWMoTWVhbk5EVkkpKSAlPiUgCiAgZ3JvdXBfYnkoRGVsdGEpICU+JSAKICBzbGljZSh3aGljaC5tYXgoTWVhbk5EU1NJKSkgJT4lIAogIHJlbmFtZShNYXhNZWFuTkRTU0ltb250aCA9IG1vbnRoLCBNYXhNZWFuTkRTU0kgPSBNZWFuTkRTU0kpCgpEZWx0YU1pbk5EVkkgPC0gCiAgRGVsdGFNZWFucyAlPiUgCiAgZmlsdGVyKHN1cmZhY2UgPT0gJ0xhbmQnKSAgJT4lIAogIHNlbGVjdCAoLWMoTWVhbk5EU1NJKSkgJT4lIAogIGdyb3VwX2J5KERlbHRhKSAlPiUgCiAgc2xpY2Uod2hpY2gubWluKE1lYW5ORFZJKSkgJT4lIAogIHJlbmFtZShNaW5NZWFuTkRWSW1vbnRoID0gbW9udGgsIE1pbk1lYW5ORFZJID0gTWVhbk5EVkkpCgpEZWx0YU1pbk5EU1NJIDwtIAogIERlbHRhTWVhbnMgJT4lIAogIGZpbHRlcihzdXJmYWNlID09ICdXYXRlcicpICAlPiUgCiAgc2VsZWN0ICgtYyhNZWFuTkRWSSkpICU+JSAKICBncm91cF9ieShEZWx0YSkgJT4lIAogIHNsaWNlKHdoaWNoLm1pbihNZWFuTkRTU0kpKSAlPiUgCiAgcmVuYW1lKE1pbk1lYW5ORFNTSW1vbnRoID0gbW9udGgsIE1pbk1lYW5ORFNTSSA9IE1lYW5ORFNTSSkKCgojam9pbiBpbnRvIDEgZGF0YWZyYW1lCkRlbHRhTWF4TWluIDwtIGxlZnRfam9pbihEZWx0YU1heE5EVkksIERlbHRhTWF4TkRTU0ksIGJ5ID0gJ0RlbHRhJykgJT4lIAogIGxlZnRfam9pbiguLERlbHRhTWluTkRWSSwgYnkgPSAnRGVsdGEnKSAlPiUgCiAgbGVmdF9qb2luKC4sRGVsdGFNaW5ORFNTSSwgYnkgPSAnRGVsdGEnKQoKI3JlbW92ZSBpbnRlcm1lZGlhdGUgZGF0YQpybShEZWx0YU1heE5EVkksIERlbHRhTWF4TkRTU0ksIERlbHRhTWluTkRTU0ksRGVsdGFNaW5ORFZJKQoKYGBgCgoKQW5kIG5vdyB3ZSBjYW4gbG9vayBhdCB0aGUgcGhhc2Ugc2hpZnRzIGJldHdlZW4gdGhlc2UgdHdvIHRpbWUgc2VyaWVzICh0aGUgdGltZXNlcmllcyBvZiBtZWFuIE5EVkkgYW5kIG1lYW4gTkRTU0kpLiBIZXJlIGFyZSB0aGUgcGhhc2Ugc2hpZnRzIChpbiBtb250aCkgZm9yIGVhY2ggZGVsdGE6CmBgYHtyIH0KI2NvbXBhcmUgb2Zmc2V0CkRlbHRhTWF4TWluIDwtIG11dGF0ZShEZWx0YU1heE1pbiwgCiAgICAgICAgICAgICAgICAgICAgICBNaW5PZmZzZXQgPSBpZl9lbHNlKGFicyhNaW5NZWFuTkRWSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpID4gNiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyIC0gYWJzKE1pbk1lYW5ORFZJbW9udGggLSBNaW5NZWFuTkRTU0ltb250aCksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGFicyhNaW5NZWFuTkRWSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpKSwKICAgICAgICAgICAgICAgICAgICAgIE1heE9mZnNldCA9IGlmX2Vsc2UoYWJzKE1heE1lYW5ORFZJbW9udGggLSBNYXhNZWFuTkRTU0ltb250aCkgPiA2LCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAxMiAtIGFicyhNYXhNZWFuTkRWSW1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhYnMoTWF4TWVhbk5EVkltb250aCAtIE1heE1lYW5ORFNTSW1vbnRoKSksCiAgICAgICAgICAgICAgICAgICAgICBPZmZzZXREaWZmID0gYWJzKE1heE9mZnNldCAtIE1pbk9mZnNldCksCiAgICAgICAgICAgICAgICAgICAgICByYW5nZU5EVkkgPSAoTWF4TWVhbk5EVkkgLSBNaW5NZWFuTkRWSSksIAogICAgICAgICAgICAgICAgICAgICAgcmFuZ2VORFNTSSA9IChNYXhNZWFuTkRTU0kgLSBNaW5NZWFuTkRTU0kpLAogICAgICAgICAgICAgICAgICAgICAgaGFsZlBlcmlvZE5EVkkgPSBpZl9lbHNlKGFicyhNYXhNZWFuTkRWSW1vbnRoIC0gTWluTWVhbk5EVkltb250aCkgPiA2LCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgMTIgLSBhYnMoTWF4TWVhbk5EVkltb250aCAtIE1pbk1lYW5ORFZJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhYnMoTWF4TWVhbk5EVkltb250aCAtIE1pbk1lYW5ORFZJbW9udGgpKSwKICAgICAgICAgICAgICAgICAgICAgIGhhbGZQZXJpb2RORFNTSSA9IGlmX2Vsc2UoYWJzKE1heE1lYW5ORFNTSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpID4gNiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyIC0gYWJzKE1heE1lYW5ORFNTSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhYnMoTWF4TWVhbk5EU1NJbW9udGggLSBNaW5NZWFuTkRTU0ltb250aCkpLCApCgojIERlbHRhTWF4TWluIDwtIAojICAgRGVsdGFNYXhNaW4gICAlPiUKIyAgIHNlbGVjdCAoYyhEZWx0YSwgTWluT2Zmc2V0LCBNYXhPZmZzZXQsIE9mZnNldERpZmYsIHJhbmdlTkRWSSwgcmFuZ2VORFNTSSxNYXhNZWFuTkRTU0ksTWluTWVhbk5EU1NJLE1heE1lYW5ORFZJLE1pbk1lYW5ORFZJKSkgCgpEZWx0YU1heE1pbgoKZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IERlbHRhLCB4ID0gTWF4T2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMTo2KSwgYnJlYWtzID0gYygxOjYpKSArCiAgZXhwYW5kX2xpbWl0cyh4ID0gYygwLDYpKSAgKyAKICBnZ3RpdGxlKCJNYXhPZmZzZXQiKQoKZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IERlbHRhLCB4ID0gTWluT2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMTo2KSwgYnJlYWtzID0gYygxOjYpKSArCiAgZXhwYW5kX2xpbWl0cyh4ID0gYygwLDYpKSAgKyAKICBnZ3RpdGxlKCJNaW5PZmZzZXQiKQoKZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IERlbHRhLCB4ID0gT2Zmc2V0RGlmZikpICsgZ2VvbV9wb2ludCgpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6NiksIGJyZWFrcyA9IGMoMTo2KSkgKwogIGV4cGFuZF9saW1pdHMoeCA9IGMoMCw2KSkgICsgCiAgZ2d0aXRsZSgiT2Zmc2V0IERpZmZlcmVuY2UiKQoKYGBgCgpgYGB7cn0KZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeCA9IE1heE1lYW5ORFZJbW9udGgpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuNSkgKyAKICBzY2FsZV95X2NvbnRpbnVvdXMoTlVMTCwgYnJlYWtzID0gTlVMTCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMToxMiksIGJyZWFrcyA9IGMoMToxMikpICsgCiAgbGFicyh4ID0gIk1vbnRoIikgKwogIGdndGl0bGUoIk1vbnRoIG9mIG1heGltdW0gbWVhbiBORFZJIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBNYXhNZWFuTkRTU0ltb250aCkpICsgCiAgZ2VvbV9kb3RwbG90KGJpbndpZHRoID0gMSxkb3RzaXplID0gMC41KSArIAogIHNjYWxlX3lfY29udGludW91cyhOVUxMLCBicmVha3MgPSBOVUxMKSArIAogIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjEyKSwgYnJlYWtzID0gYygxOjEyKSkgKyAKICBsYWJzKHggPSAiTW9udGgiKSArCiAgZ2d0aXRsZSgiTW9udGggb2YgbWF4aW11bSBtZWFuIE5EU1NJIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBNYXhPZmZzZXQpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoIk1vbnRocyBPZmZzZXQgYmV0d2VlbiBORFZJIGFuZCBORFNTSSIpCgpnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh4ID0gaGFsZlBlcmlvZE5EVkkpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoImhhbGYgcGVyaW9kIE5EVkkgIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBoYWxmUGVyaW9kTkRTU0kpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoImhhbGYgcGVyaW9kIE5EU1NJIikKCmBgYAoKYGBge3J9CiNleHRyYWN0IE5EVkkgdmFsdWUgZm9yIGVhY2ggZGVsdGEgYSB0aGUgbW9udGggb2YgbWF4IE5EU1NJIHZhbHVlCgpEZWx0YU5EVklhdE1heE5EU1NJIDwtIERlbHRhTWF4TWluICU+JQogIHNlbGVjdChEZWx0YSxNYXhNZWFuTkRTU0ltb250aCkKCkRlbHRhTWVhbnNUb0pvaW4gPC0gRGVsdGFNZWFucyAlPiUKICBmaWx0ZXIoc3VyZmFjZSA9PSAnTGFuZCcpCgpEZWx0YU5EVklhdE1heE5EU1NJIDwtIGxlZnRfam9pbihEZWx0YU5EVklhdE1heE5EU1NJLCBEZWx0YU1lYW5zVG9Kb2luLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYnkgPSBjKCdEZWx0YScsICdNYXhNZWFuTkRTU0ltb250aCcgPSdtb250aCcpKQogCkRlbHRhTkRWSWF0TWF4TkRTU0kgPC0gRGVsdGFORFZJYXRNYXhORFNTSSAlPiUKICBzZWxlY3QgKC1jKHN1cmZhY2UsIE1lYW5ORFNTSSkpCgpEZWx0YU5EVklhdE1heE5EU1NJCgpnZ3Bsb3QoRGVsdGFORFZJYXRNYXhORFNTSSwgYWVzKHggPSBNZWFuTkRWSSkpICsgCiAgZ2VvbV9kb3RwbG90KGJpbndpZHRoID0gMC4xLCBkb3RzaXplID0gMC4yNSkgKyAKICBzY2FsZV95X2NvbnRpbnVvdXMoTlVMTCwgYnJlYWtzID0gTlVMTCkgICsgCiAgeGxpbSgwLDEpICsKICBsYWJzKHggPSAiTkRWSSIpICsKICBnZ3RpdGxlKCJORFZJIGF0IG1vbnRoIG9mIG1heGltdW0gbWVhbiBORFNTSSIpCmBgYApgYGB7cn0KI2V4dHJhY3QgTkRWSSB2YWx1ZSBmb3IgZWFjaCBkZWx0YSBhIHRoZSBtb250aCBvZiBtaW4gTkRTU0kgdmFsdWUKCkRlbHRhTkRWSWF0TWluTkRTU0kgPC0gRGVsdGFNYXhNaW4gJT4lCiAgc2VsZWN0KERlbHRhLE1pbk1lYW5ORFNTSW1vbnRoKQoKRGVsdGFNZWFuc1RvSm9pbiA8LSBEZWx0YU1lYW5zICU+JQogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykKCkRlbHRhTkRWSWF0TWluTkRTU0kgPC0gbGVmdF9qb2luKERlbHRhTkRWSWF0TWluTkRTU0ksIERlbHRhTWVhbnNUb0pvaW4sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBieSA9IGMoJ0RlbHRhJywgJ01pbk1lYW5ORFNTSW1vbnRoJyA9J21vbnRoJykpCiAKRGVsdGFORFZJYXRNaW5ORFNTSSA8LSBEZWx0YU5EVklhdE1pbk5EU1NJICU+JQogIHNlbGVjdCAoLWMoc3VyZmFjZSwgTWVhbk5EU1NJKSkKCkRlbHRhTkRWSWF0TWF4TkRTU0kKCmdncGxvdChEZWx0YU5EVklhdE1pbk5EU1NJLCBhZXMoeCA9IE1lYW5ORFZJKSkgKyAKICBnZW9tX2RvdHBsb3QoYmlud2lkdGggPSAwLjEsIGRvdHNpemUgPSAwLjI1KSArIAogIHNjYWxlX3lfY29udGludW91cyhOVUxMLCBicmVha3MgPSBOVUxMKSAgKyAKICB4bGltKDAsMSkgKwogIGxhYnMoeCA9ICJORFZJIikgKwogIGdndGl0bGUoIk5EVkkgYXQgbW9udGggb2YgbWluIG1lYW4gTkRTU0kiKQpgYGAKCiAgIEFuZCBoZXJlIGFyZSB0aGUgcGhhc2Ugc2hpZnRzL29mZnNldHMgYWdhaW5zdCBvdGhlciBtZWFzdXJlZCBwYXJhbWV0ZXJzIGZvciBlYWNoIGRlbHRhLiBUaGUgcmFuZ2UsIG1heCBhbmQgbWVhbiBvZiBORFZJIGFuZCBORFNTSSBpcyBjYWxjdWxhdGVkIGZyb20gdGhlIHRpbWVzZXJpZXMsIHNvIGl0IGlzIHJlYWxseSB0aGUgbWF4LCBtaW4sIGFuZCByYW5nZSBvZiB0aGUgbW9udGhseSBtZWFucyAoaS5lLiwgdGhlIG1heGltdW0gb2YgdGhlIG1lYW5zLCB0aGUgbWluaW11bSBvZiB0aGUgbWVhbnMsIGFuZCB0aGUgcmFuZ2Ugb2YgdGhlIG1lYW4pLiBPZmZzZXQgaXMgbWVhc3VyZWQgaW4gbW9udGhzLgoKYGBge3J9CnNsaWNlMSA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gcmFuZ2VORFZJLCB4ID0gTWF4T2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgCiMgKyBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMTo2KSwgYnJlYWtzID0gYygxOjYpKSArIGV4cGFuZF9saW1pdHMoeCA9IGMoMSw2KSkgCgpzbGljZTIgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IHJhbmdlTkRWSSwgeCA9IHJhbmdlTkRTU0kpKSArIGdlb21fcG9pbnQoKSAKc2xpY2UzIDwtIGdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSByYW5nZU5EU1NJLCB4ID0gTWF4T2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgCgpzbGljZTQgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFZJLCB4ID0gcmFuZ2VORFZJKSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlNSA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gTWF4TWVhbk5EVkksIHggPSByYW5nZU5EU1NJKSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlNiA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gTWF4TWVhbk5EVkksIHggPSBNYXhPZmZzZXQpKSArIGdlb21fcG9pbnQoKSAKCnNsaWNlNyA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gTWF4TWVhbk5EU1NJLCB4ID0gTWF4TWVhbk5EVkkpKSArIGdlb21fcG9pbnQoKSAKc2xpY2U4IDwtIGdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBNYXhNZWFuTkRTU0ksIHggPSByYW5nZU5EU1NJKSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlOSA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gTWF4TWVhbk5EU1NJLCB4ID0gTWF4T2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlMTAgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFNTSSwgeCA9IHJhbmdlTkRWSSkpICsgZ2VvbV9wb2ludCgpIAoKZ3JpZC5hcnJhbmdlKHNsaWNlMSwgc2xpY2UyLCBzbGljZTMsIHNsaWNlNCwgc2xpY2U1LCBzbGljZTYsIG5jb2w9MykKZ3JpZC5hcnJhbmdlKHNsaWNlNywgc2xpY2U4LCBzbGljZTksIHNsaWNlMTAsIG5jb2w9MikKCiNyZW1vdmUgdGhvc2UgcGFuZWxzCnJtKHNsaWNlMSwgc2xpY2UyLCBzbGljZTMsIHNsaWNlNCwgc2xpY2U1LCBzbGljZTYsc2xpY2U3LCBzbGljZTgsIHNsaWNlOSwgc2xpY2UxMCkKCmBgYAoKCkpvaW4gTGF0aXR1ZGUgYW5kIGxvbmdpdHVkZSBkYXRhCgpgYGB7cn0KRGVsdGFEYXRhd0xvY2F0aW9ucyA8LSBsZWZ0X2pvaW4oRGVsdGFNYXhNaW4sIERlbHRhTG9jYXRpb25zLCBieSA9IGMoIkRlbHRhIiA9ICJEZWx0YXMiKSkKCkRlbHRhRGF0YXdMb2NhdGlvbnMgPC0gRGVsdGFEYXRhd0xvY2F0aW9ucyAlPiUKICBtdXRhdGUoQWJzb2x1dGVfTGF0aXR1ZGU9IGFicyhMYXQpKQpgYGAKCnBsb3QgcGFyYW1zIHZzIGxhdApgYGB7cn0KbG9jMSA8LSBnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHkgPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeCA9IE1heE9mZnNldCkpICsgZ2VvbV9wb2ludCgpIApsb2MyIDwtIGdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeSA9IEFic29sdXRlX0xhdGl0dWRlLCB4ID0gcmFuZ2VORFNTSSkpICsgZ2VvbV9wb2ludCgpIApsb2MzIDwtIGdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeSA9IEFic29sdXRlX0xhdGl0dWRlLCB4ID0gcmFuZ2VORFZJKSkgKyBnZW9tX3BvaW50KCkgCmxvYzQgPC0gZ2dwbG90KERlbHRhRGF0YXdMb2NhdGlvbnMsIGFlcyh5ID0gQWJzb2x1dGVfTGF0aXR1ZGUsIHggPSBNYXhNZWFuTkRWSSkpICsgZ2VvbV9wb2ludCgpIApsb2M1IDwtIGdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeSA9IEFic29sdXRlX0xhdGl0dWRlLCB4ID0gTWF4TWVhbk5EU1NJKSkgKyBnZW9tX3BvaW50KCkgCgpncmlkLmFycmFuZ2UobG9jMSwgbG9jMiwgbG9jMywgbG9jNCwgbG9jNSwgbmNvbD0yKQoKbG9jMQojZ2dzYXZlKCJsb2MxLnBkZiIsIHdpZHRoID0gNCwgaGVpZ2h0ID0gNCkKbG9jMgojZ2dzYXZlKCJsb2MyLnBkZiIsIHdpZHRoID0gNCwgaGVpZ2h0ID0gNCkKbG9jMwojZ2dzYXZlKCJsb2MzLnBkZiIsIHdpZHRoID0gNCwgaGVpZ2h0ID0gNCkKCiNyZW1vdmUgdGhvc2UgcGFuZWxzCnJtKGxvYzEsIGxvYzIsIGxvYzMsIGxvYzQsIGxvYzUpCmBgYAoKYGBge3J9CiNmaW5kIHRoZSBsaW5lYXIgbW9kZWwgCkRlbHRhT2Zmc2V0X2xtIDwtIGxtKCBBYnNvbHV0ZV9MYXRpdHVkZSB+IE1heE9mZnNldCwgZGF0YSA9IERlbHRhRGF0YXdMb2NhdGlvbnMpIAoKc3VtbWFyeShEZWx0YU9mZnNldF9sbSkKCmdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeCA9IEFic29sdXRlX0xhdGl0dWRlLCB5ID0gTWF4T2Zmc2V0KSkgKyAKICBnZW9tX3BvaW50KCkgKwogIGdlb21fc21vb3RoKG1hcHBpbmcgPSBhZXMoeCA9IEFic29sdXRlX0xhdGl0dWRlLCB5ID0gTWF4T2Zmc2V0LCApLCBtZXRob2Q9bG0gKSAKCgoKYGBgCgoKTm93IGZvciB0aGUgbWFwczoKCmBgYHtyfQp3b3JsZCA8LSBnZ3Bsb3QoKSArCiAgYm9yZGVycygid29ybGQiLCBjb2xvdXIgPSAiZ3JheTg1IiwgZmlsbCA9ICJncmF5ODAiKSArCiAgdGhlbWVfbWFwKCkgCgpEZWx0YU9mZnNldE1hcCA8LSB3b3JsZCArCiAgZ2VvbV9wb2ludChhZXMoeCA9IExvbiwgeSA9IExhdCwgY29sb3IgPSBNYXhPZmZzZXQpLAogICAgICAgICAgICAgZGF0YSA9IERlbHRhRGF0YXdMb2NhdGlvbnMsIAogICAgICAgICAgICAgc2l6ZSA9IDUpICsgc2NhbGVfY29sb3JfZ3JhZGllbnQoIGhpZ2ggPSAicmVkIiwgbG93ICA9ICJ5ZWxsb3ciKSArCiAgZ2d0aXRsZSgiT2Zmc2V0IEJldHdlZW4gTkRWSSBwZWFrIG9uIExhbmQgYW5kIE5EU1NJIHBlYWsgaW4gd2F0ZXIiKQoKRGVsdGFPZmZzZXRNYXAKI2dnc2F2ZSgiRGVsdGFPZmZzZXRNYXAucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQoKYGBgCgoKYGBge3J9CndvcmxkIDwtIGdncGxvdCgpICsKICBib3JkZXJzKCJ3b3JsZCIsIGNvbG91ciA9ICJncmF5ODUiLCBmaWxsID0gImdyYXk4MCIpICsKICB0aGVtZV9tYXAoKSAKCkRlbHRhTkRWSXJhbmdlTWFwIDwtIHdvcmxkICsKICBnZW9tX3BvaW50KGFlcyh4ID0gTG9uLCB5ID0gTGF0LCBjb2xvciA9IHJhbmdlTkRWSSksCiAgICAgICAgICAgICBkYXRhID0gRGVsdGFEYXRhd0xvY2F0aW9ucywKICAgICAgICAgICAgIHNpemUgPSA1KSArIHNjYWxlX2NvbG9yX2dyYWRpZW50KCBoaWdoID0gInJlZCIsIGxvdyAgPSAieWVsbG93IikgKyAKICBnZ3RpdGxlKCJORFZJIHJhbmdlIikKCgpEZWx0YU5EU1NJcmFuZ2VNYXAgIDwtIHdvcmxkICsKICBnZW9tX3BvaW50KGFlcyh4ID0gTG9uLCB5ID0gTGF0LCBjb2xvciA9IHJhbmdlTkRTU0kpLAogICAgICAgICAgICAgZGF0YSA9IERlbHRhRGF0YXdMb2NhdGlvbnMsIAogICAgICAgICAgICAgc2l6ZSA9IDUpICsgc2NhbGVfY29sb3JfZ3JhZGllbnQoIGhpZ2ggPSAicmVkIiwgbG93ID0gInllbGxvdyIpICsgZ2d0aXRsZSgiTkRTU0kgcmFuZ2UiKSAKCkRlbHRhTkRWSXJhbmdlTWFwIAojZ2dzYXZlKCJEZWx0YU5EVklyYW5nZU1hcC5wZGYiLCB3aWR0aCA9IDYsIGhlaWdodCA9IDQpCkRlbHRhTkRTU0lyYW5nZU1hcAojZ2dzYXZlKCJEZWx0YU5EU1NJcmFuZ2VNYXAucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQoKYGBgCgpMZXQncyBsb29rIGF0IHNvbWUgb2YgdGhlIGRhdGEuIApUbyBxdWFudGlmeSB0aGUgd2F0ZXIsIHdlIHVzZSBORFNTSS4gdG8gcXVhbnRpZnkgbGFuZCwgd2UgdXNlIE5EVkkuaGVyZSBpcyB0aGUgdGltZSBzZXJpZXMgZm9yIHRoZSBQYXJhbmEuIAoKYGBge3J9CkRlbHRhUGxvdHRlciA8LSBmdW5jdGlvbihEZWx0YU5hbWUpIHsKICAjQ291bnRzIGVhY2ggbW9udGgKICBudW1WZWcgPC0gRGVsdGFzQ2xlYW5lciAlPiUKICAgIHNlbGVjdChEZWx0YSwgc3VyZmFjZSwgbW9udGgsIG5kdmkpICU+JQogICAgZmlsdGVyKERlbHRhID09IERlbHRhTmFtZSAmIHN1cmZhY2UgPT0gIkxhbmQiICYgIWlzLm5hKG5kdmkpKSAlPiUKICAgIGdyb3VwX2J5KG1vbnRoKSAlPiUKICAgIHN1bW1hcml6ZShuID0gbigpKQogIAogIG51bVNlZCA8LSBEZWx0YXNDbGVhbmVyICU+JQogICAgc2VsZWN0KERlbHRhLCBzdXJmYWNlLCBtb250aCwgbmRzc2kpICU+JQogICAgZmlsdGVyKERlbHRhID09IERlbHRhTmFtZSAmCiAgICAgICAgICAgICBzdXJmYWNlID09ICJXYXRlciIgJiAhaXMubmEobmRzc2kpKSAlPiUKICAgIGdyb3VwX2J5KG1vbnRoKSAlPiUKICAgIHN1bW1hcml6ZShuID0gbigpKQogIAogICNIaWdobGlnaHQgdGhlIE1heGltdW0gYW5kIE1pbmltdW0gTW9udGggZm9yIGVhY2ggZGVsdGEsIE5EVkkgYW5kIE5EU1NJCiAgCiAgI0xBTkQKICBWZWcgPC0KICAgIGdncGxvdChkYXRhID0gZmlsdGVyKERlbHRhc0NsZWFuZXIsIERlbHRhID09IERlbHRhTmFtZSAmCiAgICAgICAgICAgICAgICAgICAgICAgICAgIHN1cmZhY2UgPT0gIkxhbmQiKSkgKwogICAgZ2VvbV9ib3hwbG90KGFlcyh4ID0gbW9udGgsIHkgPSBuZHZpLCBncm91cCA9IG1vbnRoKSkgKwogICAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6MTIpLCBicmVha3MgPSBjKDE6MTIpKSArCiAgICBleHBhbmRfbGltaXRzKHggPSBjKDEsIDEyKSkgKwogICAgZ2d0aXRsZShEZWx0YU5hbWUpICsKICAgICNnZW9tX3RleHQoZGF0YSA9IG51bVZlZywgYWVzKHkgPSAxLjA1LCB4ID0gbW9udGgsIGxhYmVsID0gbikpICsKICAgIGdlb21fYm94cGxvdCgKICAgICAgZGF0YSA9IGZpbHRlcigKICAgICAgICBEZWx0YXNDbGVhbmVyLAogICAgICAgIERlbHRhID09IERlbHRhTmFtZSAmCiAgICAgICAgICBzdXJmYWNlID09ICJMYW5kIiAmIG1vbnRoID09IERlbHRhTWF4TWluJE1heE1lYW5ORFZJbW9udGhbRGVsdGFNYXhNaW4kRGVsdGEgPT0gRGVsdGFOYW1lXSAKICAgICAgKSwKICAgICAgYWVzKHggPSBtb250aCwgeSA9IG5kdmksIGdyb3VwID0gbW9udGgpLAogICAgICBmaWxsID0gImdyZWVuIgogICAgKSArCiAgICBnZW9tX2JveHBsb3QoCiAgICAgIGRhdGEgPSBmaWx0ZXIoCiAgICAgICAgRGVsdGFzQ2xlYW5lciwKICAgICAgICBEZWx0YSA9PSBEZWx0YU5hbWUgJiAKICAgICAgICAgIHN1cmZhY2UgPT0gIkxhbmQiICYgbW9udGggPT0gRGVsdGFNYXhNaW4kTWluTWVhbk5EVkltb250aFtEZWx0YU1heE1pbiREZWx0YSA9PURlbHRhTmFtZV0KICAgICAgKSwKICAgICAgYWVzKHggPSBtb250aCwgeSA9IG5kdmksIGdyb3VwID0gbW9udGgpLAogICAgICBmaWxsID0gImJsdWUiCiAgICApCiAgCiAgCiAgU2VkIDwtCiAgICBnZ3Bsb3QoZGF0YSA9IGZpbHRlcihEZWx0YXNDbGVhbmVyLCBEZWx0YSA9PSBEZWx0YU5hbWUgJgogICAgICAgICAgICAgICAgICAgICAgICAgICBzdXJmYWNlID09ICJXYXRlciIpKSArCiAgICBnZW9tX2JveHBsb3QoYWVzKHggPSBtb250aCwgeSA9IG5kc3NpLCBncm91cCA9IG1vbnRoKSkgKwogICAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6MTIpLCBicmVha3MgPSBjKDE6MTIpKSArCiAgICBleHBhbmRfbGltaXRzKHggPSBjKDEsIDEyKSkgKwogICAgI2dlb21fdGV4dChkYXRhID0gbnVtU2VkLCBhZXMoeSA9IDEuMDUsIHggPSBtb250aCwgbGFiZWwgPSBuKSkgKwogICAgZ2VvbV9ib3hwbG90KAogICAgICBkYXRhID0gZmlsdGVyKAogICAgICAgIERlbHRhc0NsZWFuZXIsCiAgICAgICAgRGVsdGEgPT0gRGVsdGFOYW1lICYKICAgICAgICAgIHN1cmZhY2UgPT0gIldhdGVyIiAmIG1vbnRoID09IERlbHRhTWF4TWluJE1heE1lYW5ORFNTSW1vbnRoW0RlbHRhTWF4TWluJERlbHRhID09IERlbHRhTmFtZV0KICAgICAgKSwKICAgICAgYWVzKHggPSBtb250aCwgeSA9IG5kc3NpLCBncm91cCA9IG1vbnRoKSwKICAgICAgZmlsbCA9ICJncmVlbiIKICAgICkgKwogICAgZ2VvbV9ib3hwbG90KAogICAgICBkYXRhID0gZmlsdGVyKAogICAgICAgIERlbHRhc0NsZWFuZXIsCiAgICAgICAgRGVsdGEgPT0gRGVsdGFOYW1lICYKICAgICAgICAgIHN1cmZhY2UgPT0gIldhdGVyIiAmIG1vbnRoID09IERlbHRhTWF4TWluJE1pbk1lYW5ORFNTSW1vbnRoW0RlbHRhTWF4TWluJERlbHRhID09IERlbHRhTmFtZV0KICAgICAgKSwKICAgICAgYWVzKHggPSBtb250aCwgeSA9IG5kc3NpLCBncm91cCA9IG1vbnRoKSwKICAgICAgZmlsbCA9ICJibHVlIgogICAgKQogIAogIHJldHVybihncmlkLmFycmFuZ2UoVmVnLCBTZWQsIG5yb3cgPSAyKSkKfQpgYGAKCgpIZXJlIGlzIGFuICdPZmZzZXQnIGV4YW1wbGUgZnJvbSB0aGUgUGFyYW5hOgpgYGB7cn0KRGVsdGFQbG90dGVyKCJQYXJhbmEiKQpgYGAKCgphbmQgbm93IGFuIGluLXBoYXNlIGV4YW1wbGUgZnJvbQpgYGB7cn0KRGVsdGFQbG90dGVyKCJNYWdkYWxlbmEiKQpgYGAKCgphbmQgdGhlIEVicm86CmBgYHtyfQpEZWx0YVBsb3R0ZXIoIkVicm8iKQpgYGAKCgphbmQgdGhlIFRhbmEKYGBge3J9CkRlbHRhUGxvdHRlcigiVGFuYSIpCmBgYAoKYW5kIHRoZSBUYW5hCmBgYHtyfQpEZWx0YVBsb3R0ZXIoIkdvZGF2YXJpIikKRGVsdGFQbG90dGVyKCJLcmlzaG5hIikKYGBgCgoKCk5vdyBmb3Igc29tZSB3b3JrIHdpdGggR1JEQyBkYXRhOgpgYGB7cn0KI2ltcG9ydCB0aGUgZGF0YSAobW9udGhseSBtZWFucyBmb3IgMjEgc3RhdGlvbnMpCkRlbHRhc0dSREMgIDwtIHJlYWRfY3N2KCIuLi9kYXRhL0dSRENzdGF0aW9ucy5jc3YiKQoKI2NhbGN1bGF0ZSB0aGUgbWVhbiBvZiB0aGUgbW9udGhseSBtZWFucwpEZWx0YXNHUkRDd01lYW4gPC0gRGVsdGFzR1JEQyAlPiUgCiAgICByb3d3aXNlKCkgJT4lIAogICAgbXV0YXRlKE1NRD1tZWFuKGMoSmFudWFyeSxGZWJydWFyeSxNYXJjaCxBcHJpbCwKICAgICAgICAgICAgICAgICAgICBNYXksSnVuZSxKdWx5LEF1Z3VzdCwKICAgICAgICAgICAgICAgICAgICBTZXB0ZW1iZXIsT2N0b2JlcixOb3ZlbWJlcixEZWNlbWJlcikpKQoKRGVsdGFzR1JEQ3dNZWFuIDwtIERlbHRhc0dSREN3TWVhbiAlPiUgCiAgcm93d2lzZSgpICU+JSAKICBtdXRhdGUoUmFuZ2VfRGlzY2hhcmdlID0gbWF4KGMoSmFudWFyeSxGZWJydWFyeSxNYXJjaCxBcHJpbCwKICAgICAgICAgICAgICBNYXksSnVuZSxKdWx5LEF1Z3VzdCwKICAgICAgICAgICAgICBTZXB0ZW1iZXIsT2N0b2JlcixOb3ZlbWJlcixEZWNlbWJlcikpIC0gCiAgICAgICAgICAgbWluKGMoSmFudWFyeSxGZWJydWFyeSxNYXJjaCxBcHJpbCwKICAgICAgICAgICAgICBNYXksSnVuZSxKdWx5LEF1Z3VzdCwKICAgICAgICAgICAgICBTZXB0ZW1iZXIsT2N0b2JlcixOb3ZlbWJlcixEZWNlbWJlcikpCiAgICAgICAgICAgKQogICAgICAgIAoKI2pvaW4gdG8gbG9jYXRpb24gZGF0YToKRGVsdGF3TG9jR1JEQyA8LSBsZWZ0X2pvaW4oRGVsdGFEYXRhd0xvY2F0aW9ucywgRGVsdGFzR1JEQ3dNZWFuLCBieSA9IGMoIkRlbHRhIiA9ICJEZWx0YXMiKSkKCgojcGxvdCBtZWFuIG9mIG1vbnRobHkgbWVhbnMgYWdhaW5zdCBORFNTSQpnZ3Bsb3QoRGVsdGF3TG9jR1JEQywgYWVzKHkgPSBSYW5nZV9EaXNjaGFyZ2UsIHggPSByYW5nZU5EU1NJKSkgKyBnZW9tX3BvaW50KCkgKyBzY2FsZV95X2NvbnRpbnVvdXModHJhbnM9J2xvZzEwJykKCiNnZ3NhdmUoIkdSRENORFNTSS5wZGYiLCB3aWR0aCA9IDYsIGhlaWdodCA9IDQpCgpnZ3Bsb3QoRGVsdGF3TG9jR1JEQywgYWVzKHkgPSBSYW5nZV9EaXNjaGFyZ2UsIHggPSBNYXhNZWFuTkRTU0kpKSArIGdlb21fcG9pbnQoKSArIHNjYWxlX3lfY29udGludW91cyh0cmFucz0nbG9nMTAnKQpgYGAKYGBge3J9CiNpbXBvcnQgdGhlIGRhdGEgKG1vbnRobHkgbWVhbnMgZm9yIDIxIHN0YXRpb25zKQpEZWx0YXNHUkRDICA8LSByZWFkX2NzdigiLi4vZGF0YS9HUkRDc3RhdGlvbnMuY3N2IikKYGBgCgoKYGBge3J9CiNjYWxjdWxhdGUgdGhlIG1lYW4gb2YgdGhlIG1vbnRobHkgbWVhbnMKRGVsdGFzR1JEQ3dNZWFuIDwtIERlbHRhc0dSREMgJT4lIAogICAgcm93d2lzZSgpICU+JSAKICAgIG11dGF0ZShNTUQ9bWVhbihjKEphbnVhcnksRmVicnVhcnksTWFyY2gsQXByaWwsCiAgICAgICAgICAgICAgICAgICAgTWF5LEp1bmUsSnVseSxBdWd1c3QsCiAgICAgICAgICAgICAgICAgICAgU2VwdGVtYmVyLE9jdG9iZXIsTm92ZW1iZXIsRGVjZW1iZXIpKSkKCkRlbHRhc0dSREN3TWVhbiA8LSBEZWx0YXNHUkRDd01lYW4gJT4lIAogIHJvd3dpc2UoKSAlPiUgCiAgbXV0YXRlKFJhbmdlX0Rpc2NoYXJnZSA9IG1heChjKEphbnVhcnksRmVicnVhcnksTWFyY2gsQXByaWwsCiAgICAgICAgICAgICAgTWF5LEp1bmUsSnVseSxBdWd1c3QsCiAgICAgICAgICAgICAgU2VwdGVtYmVyLE9jdG9iZXIsTm92ZW1iZXIsRGVjZW1iZXIpKSAtIAogICAgICAgICAgIG1pbihjKEphbnVhcnksRmVicnVhcnksTWFyY2gsQXByaWwsCiAgICAgICAgICAgICAgTWF5LEp1bmUsSnVseSxBdWd1c3QsCiAgICAgICAgICAgICAgU2VwdGVtYmVyLE9jdG9iZXIsTm92ZW1iZXIsRGVjZW1iZXIpKQogICAgICAgICAgICkKICAgICAgICAKCiNqb2luIHRvIGxvY2F0aW9uIGRhdGE6CkRlbHRhd0xvY0dSREMgPC0gbGVmdF9qb2luKERlbHRhRGF0YXdMb2NhdGlvbnMsIERlbHRhc0dSREN3TWVhbiwgYnkgPSBjKCJEZWx0YSIgPSAiRGVsdGFzIikpCgoKI3Bsb3QgbWVhbiBvZiBtb250aGx5IG1lYW5zIGFnYWluc3QgTkRTU0kKZ2dwbG90KERlbHRhd0xvY0dSREMsIGFlcyh5ID0gUmFuZ2VfRGlzY2hhcmdlLCB4ID0gcmFuZ2VORFNTSSkpICsgZ2VvbV9wb2ludCgpICsgc2NhbGVfeV9jb250aW51b3VzKHRyYW5zPSdsb2cxMCcpCgojZ2dzYXZlKCJHUkRDTkRTU0kucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQpgYGAKCmBgYHtyfQpnZ3Bsb3QoRGVsdGF3TG9jR1JEQywgYWVzKHkgPSBSYW5nZV9EaXNjaGFyZ2UsIHggPSBNYXhNZWFuTkRTU0kpKSArIGdlb21fcG9pbnQoKSArIHNjYWxlX3lfY29udGludW91cyh0cmFucz0nbG9nMTAnKQpgYGAKCmBgYHtyfQojcmVuYW1lIHRoZSBtb250aHMgYnkgbnVtYmVycyBhbmQgdGlkeSB0aGUgR1JEQyBkYXRhCkRlbHRhc0Rpc2NoYXJnZSA8LSBEZWx0YXNHUkRDICU+JQogIHJlbmFtZShEZWx0YSA9IERlbHRhcywiMSIgPSBKYW51YXJ5LCAiMiI9IEZlYnJ1YXJ5LCAiMyI9IE1hcmNoLCAiNCI9IEFwcmlsLAogICAgICAgICAiNSI9TWF5LCAiNiI9SnVuZSwgIjciPUp1bHksICI4Ij0gQXVndXN0LCAiOSIgPSBTZXB0ZW1iZXIsICIxMCI9T2N0b2JlciwgCiAgICAgICAgICIxMSI9Tm92ZW1iZXIsICIxMiI9RGVjZW1iZXIpICU+JQogIHNlbGVjdChEZWx0YSwgIjEiICwgIjIiICwgIjMiLCAiNCIsIjUiLCAiNiIsICI3IiwgIjgiLCAiOSIsICIxMCIsICIxMSIsICIxMiIpICU+JQogIHBpdm90X2xvbmdlcigtRGVsdGEsIG5hbWVzX3RvID0gIm1vbnRoIiwgdmFsdWVzX3RvID0gImRpc2NoYXJnZSIpCgojZmluZCBtYXggR1JEQyBtb250aCBmb3IgZWFjaCBkZWx0YQpEZWx0YU1heERpc2NoYXJnZSA8LSAKICBEZWx0YXNEaXNjaGFyZ2UgJT4lIAogIGdyb3VwX2J5KERlbHRhKSAlPiUgCiAgc2xpY2Uod2hpY2gubWF4KGRpc2NoYXJnZSkpICU+JSAKICByZW5hbWUoTWF4RGlzY2hhcmdlTW9udGggPSBtb250aCwgTWF4RGlzY2hhcmdlID0gZGlzY2hhcmdlKSAKCgpEZWx0YU1heERpc2NoYXJnZSRNYXhEaXNjaGFyZ2VNb250aCA9IGFzLm51bWVyaWMoRGVsdGFNYXhEaXNjaGFyZ2UkTWF4RGlzY2hhcmdlTW9udGgpCgoKI2pvaW4gd2l0aCBvdGhlciBkZWx0YSBkYXRhCkRlbHRhTWF4TWluRGlzY2hhcmdlIDwtIGxlZnRfam9pbihEZWx0YU1heE1pbiwgRGVsdGFNYXhEaXNjaGFyZ2UsIGJ5ID0gJ0RlbHRhJykKCiNjYWxjdWxhdGUgb2Zmc2V0CkRlbHRhTWF4TWluRGlzY2hhcmdlIDwtIERlbHRhTWF4TWluRGlzY2hhcmdlICU+JQogIG11dGF0ZSggRGlzc09mZiA9IGlmX2Vsc2UoYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpID4gNiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyIC0gYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpKQogICAgICAgICAgKQoKI0NvbXBhcmUgb2Zmc2V0IHdpdGggTkRTU0kgKGRlbHRhbWF4bWluJG1heG1lYW5ORFNTSW1vbnRoKQpnZ3Bsb3QoRGVsdGFNYXhNaW5EaXNjaGFyZ2UsIGFlcyh5ID0gRGVsdGEsIHggPSBEaXNzT2ZmKSkgKyBnZW9tX3BvaW50KCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMTo2KSwgYnJlYWtzID0gYygxOjYpKSArCiAgZXhwYW5kX2xpbWl0cyh4ID0gYygwLDYpKSAgKyAKICBnZ3RpdGxlKCJEaXNPZmZzZXQiKQpgYGAKCmBgYHtyfQojam9pbiBsYXQgZGF0YQpEZWx0YURhdGF3TG9jYXRpb25zIDwtIGxlZnRfam9pbihEZWx0YU1heE1pbkRpc2NoYXJnZSwgRGVsdGFMb2NhdGlvbnMsIGJ5ID0gYygiRGVsdGEiID0gIkRlbHRhcyIpKQoKRGVsdGFEYXRhd0xvY2F0aW9ucyA8LSBEZWx0YURhdGF3TG9jYXRpb25zICU+JQogIG11dGF0ZShBYnNvbHV0ZV9MYXRpdHVkZT0gYWJzKExhdCkpCgojIHBsb3Qgb2Zmc2V0IG9uIGdyYXBoIGJ5IGxhdApnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHggPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeSA9IERpc3NPZmYpKSArIGdlb21fcG9pbnQoKSArCiAgc2NhbGVfY29sb3JfZ3JhZGllbnQobG93ID0gInllbGxvdyIsIGhpZ2ggPSAicmVkIiwgbmEudmFsdWUgPSBOQSkgCiAgIysgZ2VvbV9zbW9vdGgobWFwcGluZyA9IGFlcyh4ID0gQWJzb2x1dGVfTGF0aXR1ZGUsIHkgPSBEaXNzT2ZmLCApLCBtZXRob2Q9bG0gKSAKCiNnZ3NhdmUoIkRpc09mZnNldC5wZGYiLCB3aWR0aCA9IDYsIGhlaWdodCA9IDQpCmBgYAoKYGBge3J9CiNwbG90IG9mZnNldCBvbiBtYXAKRGVsdGFEaXNPZmZzZXRNYXAgPC0gd29ybGQgKwogIGdlb21fcG9pbnQoYWVzKHggPSBMb24sIHkgPSBMYXQsIGNvbG9yID0gRGlzc09mZiksCiAgICAgICAgICAgICBkYXRhID0gRGVsdGFEYXRhd0xvY2F0aW9ucywgCiAgICAgICAgICAgICBzaXplID0gNSkgKyBzY2FsZV9jb2xvcl9ncmFkaWVudCggaGlnaCA9ICJyZWQiLCBsb3cgID0gInllbGxvdyIpICsKICBnZ3RpdGxlKCJPZmZzZXQgQmV0d2VlbiBHUkRDIGRpc2NoYXJnZSBwZWFrIGFuZCBORFNTSSBwZWFrIGluIHdhdGVyIikKCkRlbHRhRGlzT2Zmc2V0TWFwCiNnZ3NhdmUoIkRlbHRhT2Zmc2V0TWFwLnBkZiIsIHdpZHRoID0gNiwgaGVpZ2h0ID0gNCkKCmBgYAoKCgo=